##########################################################################################
#### 20231211 ArchR分析步骤 ####

options(stringsAsFactors=F)
library(Seurat)
library(ArchR)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--rna_data_file"), type = "character"),
    make_option(c("--input_dir"), type = "character"),
    make_option(c("--cluster_cell_file"), type = "character"),
    make_option(c("--cpu"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
  input_dir <- "~/20231121_singleMuti/input"
  out_path <- "~/20231121_singleMuti/results/atac_res"
  rna_data_file <- "~/20231121_singleMuti/input/testis_combined.Rdata"
  cluster_cell_file <- "~/20231121_singleMuti/config/cluster_celltype.csv"
  cpu <- 10

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_data_file <- opt$rna_data_file
input_dir <- opt$input_dir
cluster_cell_file <- opt$cluster_cell_file
cpu <- as.numeric(opt$cpu)
out_path <- opt$out_path

##########################################################################################

dir.create( out_path , recursive = T )
setwd( out_path )

addArchRThreads(threads = cpu)
addArchRGenome("hg38") 

##########################################################################################
## 转录组文件
rna_data <- load(rna_data_file)

## atac文件
inputFiles <- c(
  paste0( input_dir , "/testis01_atac_fragments.tsv.gz" ) ,
  paste0( input_dir , "/testis02_atac_fragments.tsv.gz" ) ,
  paste0( input_dir , "/testis03_atac_fragments.tsv.gz" ) 
  )

## 细胞注释文件
cluster_cell <- fread(cluster_cell_file , header = F)

##########################################################################################
## 读atac的文件
names(inputFiles) <- c("testis01","testis02","testis03") ## 3个样本的fragments文件
ArrowFile <- createArrowFiles(inputFiles = inputFiles,sampleNames = names(inputFiles),minTSS = 4,minFrags = 1000,addTileMat = TRUE,addGeneScoreMat = TRUE) ## 质控
testis <- ArchRProject(ArrowFiles = ArrowFile, outputDirectory = out_path, copyArrows = TRUE)  ## outputDirectory就是用来放archr中间文件的目录

##########################################################################################
# 和表达数据Rdata中的细胞做交集,选择scRNA-seq和scATAC-seq数据都通过质控的细胞
# 例(这一步是因为如果要做addGeneIntegrationMatrix这一步的时候,细胞必须要是一样的)
scrnat <- testis_combined
scrnat$cell_type <- ifelse(scrnat$seurat_clusters==0,"Myoid cells",
  ifelse(scrnat$seurat_clusters==1,"Round&ElongateS.tids",
  ifelse(scrnat$seurat_clusters==2,"Sperm",
  ifelse(scrnat$seurat_clusters==3,"Leydig cells",
  ifelse(scrnat$seurat_clusters==4,"Differenting&Differented SPG",
  ifelse(scrnat$seurat_clusters==5,"SSC",
  ifelse(scrnat$seurat_clusters==6,"Endothelial cells",
  ifelse(scrnat$seurat_clusters==7,"Early stage of spermatids",
  ifelse(scrnat$seurat_clusters==8,"Leptotene",
  ifelse(scrnat$seurat_clusters==9,"Zygotene",
  ifelse(scrnat$seurat_clusters==10,"Diplotene",
  ifelse(scrnat$seurat_clusters==11,"Sertoli cells",
  ifelse(scrnat$seurat_clusters==12,"Patchytene",
  ifelse(scrnat$seurat_clusters==13,"Patchytene",
  ifelse(scrnat$seurat_clusters==14,"Macrophages",
  ifelse(scrnat$seurat_clusters==15,"Pericytes",
  ifelse(scrnat$seurat_clusters==16,"NKT cells","unkown")))))))))))))))))

cell_id <- data.frame(bc=colnames(scrnat))
cell_id$id2 <- gsub("_" , "#" , cell_id$bc)
cell_id$cell_type <- scrnat$cell_type
cell_id$seurat_clusters <- scrnat$seurat_clusters
rownames(cell_id) <- cell_id$id2

cellname_index <- which(rownames(testis@cellColData) %in% cell_id$id2)  

testis2c <- testis
testis2c@cellColData <- testis2c@cellColData[cellname_index,]

## RNA总共24735细胞
## atac总共47725细胞
## 交集是22261个细胞

##########################################################################################
# scATAC-seq数据的聚类
testis2c <- addIterativeLSI(ArchRProj = testis2c,useMatrix = "TileMatrix", name = "IterativeLSI", iterations = 2, clusterParams = list(resolution = c(0.2), sampleCells = 10000, n.start = 10),varFeatures = 25000, dimsToUse = 1:30)
testis2c <- addClusters(input = testis2c,reducedDims = "IterativeLSI",method = "Seurat",name = "Clusters",resolution = 0.8)
testis2c <- addUMAP(testis2c,saveModel = FALSE)

## 依据rna的数据定义细胞类型
testis2c@cellColData$cell_type <- cell_id[rownames(testis2c@cellColData),"cell_type"]
testis2c@cellColData$seurat_clusters <- cell_id[rownames(testis2c@cellColData),"seurat_clusters"]

# 调用macs2 
pathToMacs2 <- findMacs2()
#pathToMacs2 <- "~/miniconda3/envs/archr/bin/macs2"

# call peak
testis2c <- addGroupCoverages(testis2c,groupBy = "cell_type")  ## 创建伪批量重复
testis2c <- addReproduciblePeakSet(testis2c,groupBy = "cell_type")  ## 迭代重叠峰
testis2c <- addPeakMatrix(testis2c)

## 输出
testis_combined_peak <- testis2c
out_file <- paste0( out_path , "/testis_combined_peak.Rdata" ) 
save( testis_combined_peak, file = out_file )

## peak的counts矩阵输出为seruat对象
all_peak <- data.frame(testis_combined_peak@peakSet)
all_peak$cell_type <- names(testis_combined_peak@peakSet)

cell_id$seurat_clusters <- as.character(cell_id$seurat_clusters)
cell_id$seurat_clusters <- ifelse( cell_id$seurat_clusters == 12 | cell_id$seurat_clusters == 13 , 17 , cell_id$seurat_clusters  )
cell_id$seurat_clusters <- paste0( cell_id$seurat_clusters )
all_peak <- merge( all_peak , unique(cell_id[ , c("cell_type" , "seurat_clusters")]) , by = "cell_type")
all_peak <- all_peak[, c(2:ncol(all_peak) , 1) ]

out_file <- paste0( out_path , "/testis_combined_peak.tsv" ) 
write.table( all_peak , out_file , row.names = F , quote = F , sep = "\t" )

##########################################################################################
## 输出符合质控的细胞列表
cell_qc <- data.frame(V1 = rownames(testis_combined_peak@cellColData))

out_file <- paste0( out_path , "/testis_merge_all_qc_barcode_new.txt" ) 
write.table( cell_qc , out_file , row.names = F , quote = F , sep = "\t" )
